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Abstract 

A novel and efficient approach which is based on the framework of isogeometric 
analysis for elliptic homogenization problems is proposed. These problems possess highly 
oscillating coefficients leading to extremely high computational expenses while using 
traditional finite element methods. The isogeometric analysis heterogeneous multiscale 
method (IGA-HMM) investigated in this paper is regarded as an alternative approach 
to the standard Finite Element Heterogeneous Multiscale Method (FE-HMM) which 
is currently an effective framework to solve these problems. The method utilizes non- 
uniform rational B-splines (NURBS) in both macro and micro levels instead of standard 
Lagrange basis. Beside the ability to describe exactly the geometry, it tremendously 
facilitates high-order macroscopic/microscopic discretizations thanks to the flexibility of 
refinement and degree elevation with an arbitrary continuity level provided by NURBS 
basis functions. A priori error estimates of the discretization error coming from macro 
and micro meshes and optimal micro refinement strategies for macro/micro NURBS 
basis functions of arbitrary orders are derived. Numerical results show the excellent 
performance of the proposed method. 

Keywords: Isogeometric analysis (IGA), NURBS, IGA-HMM, Heterogeneous 
Multiscale Method, Homogenization. 



1. Introduction 

Homogenization is a branch of mathematics and engineering which studies partial 
differential equations (PDEs) with rapidly oscillating coefficients. These kinds of equa- 
tions describe various processes in inhomogeneous materials with rapidly oscillating mi- 
cro structures, such as composite and perforated materials, and thus play an important 
role in physics, engineering and modern technologies. The aim of homogenization is to 
"average out" the heterogeneities at micro-scale and describe the effective properties at 
macro-scale of such phenomena. In other words, we want to know the behaviors of the 
systems as homogeneous ones. In analytical approaches, such as in [7|[T3]. homogenized 
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equations are derived. However, the coefficients of these equations are only computed 
explicitly in some special cases, such as when the medium follows some periodic as- 
sumptions, and not explicitly available in general. Furthermore, full computations with 
complex scale interactions of the heterogeneous system are very ineffective due to high 
computational cost. Thus, to solve these problems, advanced computational technologies 
have been developed. 

Literature reviews on various multiscale approaches can be found in pjU UHl [14] . 
In this paper, we focus on the heterogeneous multiscale method (HMM), which was 
proposed in [20] . Reviews on HMM are presented in [2U III El US] . This method provides 
a general framework which allows ones to develop various approaches to homogenization 
problems. The simplest one is the finite element heterogeneous multiscale method (FE- 
HMM), which uses standard finite elements such as simplicial or quadrilateral ones in 
both macroscopic and microscopic level. Solving the so-called micro problems (with 
a suitable set up) in sampling domains around traditional Gauss integration points at 
macro level allows one to approximate the missing effective information for the macro 
solver. 

The standard finite element method (FEM), albeit very popular in various fields of 
sciences and engineering, still has some shortcomings which affect the efficiency of the 
FE-HMM. Firstly, the discretized geometry through mesh generation is required. This 
process often results in geometrical errors even with the higher-order FEM. Also, the 
communication between the geometry model and the mesh program during the analysis 
process is always needed and this constitutes the large part of the overall process [TT] . 
especially for industrial problems. Secondly, lower-order formulations, such as FE-HMM 
based on the four-node quadrilateral element (Q4), often require extremely fine meshes 
to produce approximate solutions with a desired accuracy for complicated problems. 
This prevents multiscale analyses from being run on personal computers. Thirdly, high- 
order formulations still put some restrictions on the element topologies (for example, 
the connection of different type of corner, center, or internal nodes) and only possess 
C° continuity. These disadvantages lead to an increase in the number of micro coupling 
problems and thus increase the overall computational cost. Hence, there is a need to 
consider alternative methods to tackle these issues. 

Among advanced numerical methods, the so-called isogeometric analysis (IGA), 
where NURBS are used as basis functions, emerges as the most potential candidate. 
The isogeometric analysis was first proposed by Hughes and co-workers [H] and now 
has attracted the attention of academic as well as engineering community all over the 
world. The IGA provides a framework in which the gap between Computer Aided De- 
sign (CAD) and Finite Element Analysis (FEA) may be reduced. This is achieved in 
IGA by employing the same basis functions to describe both the geometry of the domain 
of interest and the field variables. While the standard FEM uses basis functions which 
are based on Lagrange polynomials, isogeometric approach utilizes more general basis 
functions such as B-splines and NURBS that are commonly used in CAD geometry. 
The exact geometry is therefore maintained at the coarsest level of discretization and 
re-meshing is seamlessly performed on this coarsest level without any further communi- 
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cation with CAD geometry. Furthermore, B-splines (or NURBS) provide a flexible way 
to make refinement, de-refinement, and degree elevation [TT] . They allow us to achieve 
easily the smoothness of arbitrary continuity in comparison with the C° continuity pro- 
vided by the traditional FEM. For a reference on IGA, we recommend the excellent book 
[8] and we refer to the NURBS book [16] for a geometric description. 

In this work, we introduce a new approach: a so-called isogeometric analysis heteroge- 
neous multiscale method (IGA-HMM) which utilizes NURBS as basis functions for both 
exact geometric representation and analysis. The NURBS are used as basis functions for 
both macro and micro element spaces, where the FE-HMM only employs standard FEM 
basis. This tremendously facilitate high-order macroscopic/microscopic discretizations 
by a flexibility of refinements and degree elevations with an arbitrary continuity of basis 
functions. As will be demonstrated later in the numerical examples section, elliptic ho- 
mogenization problems can be solved, by using the proposed IGA-HMM, effectively on 
a personal computer which is not the case for FE-HMM with (bi) linear basis functions if 
high accuracies are needed. The FE-HMM with bi (linear) elements often requires very 
fine macro meshes (thus a high number of micro problems) that are only supported by 
powerful computers. It is obvious that one can use high order Lagrange elements in 
the framework of FE-HMM to achieve the same goal. However, high order Lagrange 
elements cannot be constructed as straightforwardly as NURBS and more importantly 
they are only C° continuity whereas NURBS elements are C V ~ Y (p is the NURBS order) 
by definition. Therefore, IGA-HMM is able to solve high order homogenization problems 
such as plate/shell homogenization problems. One can consider our IGA-HMM as an 
efficient high order NURBS-based FE-HMM. We refer to [T7j for a related work on high 
order FE-HMM which, however, limits to cubic macro elements and quadratic micro 
elements. A priori error estimates of the discretization error coming from macro and 
micro meshes are provided. Optimal micro refinement strategies for macro and micro 
NURBS basis functions of arbitrary orders are presented and thus an optimal value for 
the micro NURBS order is obtained as a function of the macro NURBS basis order. 
We have also observed that C°, not C p_1 , high order NURBS should be used at the 
microscale. It should be mentioned that the proposed IGA-HMM is very similar to the 
standard FE-HMM, this means that the implementation is very simple, any existing 
coding framework of FE-HMM e.g. the one given in [4] can be readily re-used. 

The paper is arranged as follows: a brief introduction to the B-splines, NURBS and 
IGA is given in Section 2. Section 3 describes the isogeometric analysis heterogeneous 
multiscale method. Numerical examples are provided in Section 4 and Section 5 closes 
the paper with some concluding remarks. 

2. NURBS-based isogeometric analysis fundamentals 

2.1. Knot vectors and basis functions 

Univariate B-spline basis functions are constructed from knot vectors. Given two 
positive integers p and n and let 5 = £2, £ n +p+i] be a non-decreasing sequence 
of parameter values, & < = 1, + p. The & are called knots, S is the set of 
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coordinates in the parametric space, p is the polynomial degree, and n is the number 
of basis functions. If all knots are equally spaced the knot vector is called uniform. 
Otherwise, they are called a non-uniform knot vector. When the first and the last knots 
are repeated p+1 times, the knot vector is called open. An important property of open 
knot vectors is that the resulting basis functions are interpolatory at the ends of the 
parametric space. A B-spline basis function is C°° continuous inside knot spans and 
C p ~ l continuous at a single knot. If an interior knot value repeats one more time, it is 
then called a multiple knot. At a knot of multiplicity k the continuity is C p ~ k . Given 
a knot vector, the B-spline basis functions A^ p (£) are defined starting with the zeroth 
order (p = 0) basis functions 



Ni,o(0 = 

and for p > 1, they are defined recursively as follows [16 
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For p = and p = 1, the NURBS basis functions of isogeometric analysis are identical 
to those of standard piecewise constant and linear finite elements, respectively. Never- 
theless, for p > 2, they are different [11]. Therefore, the present work will consider the 
basis functions with p > 2. 

Fig. [T] illustrates a set of univariate quadratic, cubic and quartic B-spline basis func- 
tions for open uniform knot vectors S = {0, 0, 0, |, 1, 1, 1}, 5 = {0, 0, 0, 0, |, 1, 1, 1, 1} 
and 5 = {0, 0, 0, 0, 0, |, 1, 1, 1, 1, 1}, respectively. 






Figure 1: An illustration of quadratic, cubic and quartic B-spline basis functions. 



2.2. NURBS curves and surfaces 
The B-spline curve is defined as 

n 

c(o=x;^(OPi, (3) 

1=1 

where are the control points and N^ v (£) is the p^-degree B-spline basis function 
defined on an open knot vector. 
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The B-spline surfaces are defined by the tensor product of the basis functions in 
two parametric dimensions £ and 77 with two knot vectors S = {£1, £2..., £n+p+i} an d 
ft = {??i,%...,Wg+l} as 

n m 

s(^) = EE ^ (o M « fa) p « > ( 4 ) 

i=i j=i 

where is the bidirectional control net, iV^ (£) and M^ q {rj) are the B-spline basis 
functions defined on the knot vectors over an n x m net of control points T?ij. Following 
the traditional FEM notation, Eq. Q is rewritten as follows 

nxm 

S(^ V ) = Y / ^(^v)Pi, (5) 

where Nj (£,77) = Ni^ p (£) M^ q (rj) is the bivariate B-spline basis function associated with 
node I. Fig. [2] plots the bivariate quadratic and cubic B-spline basis functions. 




Figure 2: Bivariate quadratic and cubic B-spline basis functions with open uniform knot 
vectors 5 = U = {0, 0, 0, §, 1, 1, 1}, S = U = {0, 0, 0, 0, §, 1, 1, 1, 1}. 

B-splines are convenient for free-form modelling, but they lack the ability to ex- 
actly represent some simple shapes like circles and ellipsoids. Non-uniform rational 
B-splines (NURBS) extend B-splines since they allow exact representation of conic sec- 
tions. NURBS are obtained by augmenting every point in the control mesh Pj with the 
weight wj. The weighting function is constructed as follows 

nxm 

v>(t,v) = J2Ni(t,v)v>i- (6) 

1=1 

The NURBS surfaces are then defined by 

nxm 

E N» (£, V ) 

^/Pj nxm 
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where Nj (£, rj) = ^r^p are the rational basis functions. An example of quadratic 
NURBS curve and surface is demonstrated in Fig. [3j It is clear that, due to the par- 
tition of unity property of B-splines basis, if all the weights are units then NURBS 
curves/surfaces become B-splines curves /surfaces. 




Figure 3: A quadratic NURBS curve (left) and a quadratic NURBS surface with physical 
mesh and control grid (right) 



2.3. Refinement 

In addition to the classical /i-refinement (knot insertion in CAD notation) and p- 
refinement (degree elevation in CAD notation) supported by FEM, IGA has a new, 
unique, and more economical version: fc-refinement which is a process in which degree 
elevation is followed by knot insertion. The fc-refinement allows us to elevate the degree 
of NURBS objects while also obtaining higher continuities across element boundaries, 
whereas the standard FEM always stay at C°. In comparison to p- version, k- version 
require much less number of basis functions while maintaining the same convergence 
rate. The readers are referred to [HI [8] for more details. 

2.4- Isoparametric concept 

Isogeometric analysis also employs the isoparametric concept as in standard FEM- 
the same NURBS basis functions are used for both the description of the geometry and 
the unknown field: 

nxm nxm 

x(£, V ) = Ni (£, v)Pi, u H (x($, 77)) = J2 N i & V)ui, (8) 

where nxm represent the number of basis functions, x T = (x, y) is the physical coordi- 
nates vector, Ni (£, rj) is the NURBS basis function and ui is the degree of freedom (dof) 
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of u H (which can be the displacements for solid mechanics problems, the temperatures 
for thermal problems etc.) at the control point /, respectively. The control points play 
the role of nodes in a FEM setting and the IGA elements are defined as the non-zero 
knot spans. 

2.5. Properties of IGA 

The NURBS-based IGA and FEA share a lot of commons such as [llj: compact 
support, partition of unity, affine covariance, isoparametric concept, patch test satisfied. 
These properties allow implementations in the IGA to follow the framework of Galerkin 
method as in the FEA. On the other hand, there are important properties which IGA 
possesses but FEA does not. Some differences between NURBS-based IGA and FEA 
are presented in Tab. [TJ Note that one weak point of NURBS-based IGA is the lack of 
interpolation property which renders the application of Dirichlet boundary conditions a 
bit more involved than traditional FEM. We refer to pj] for details on this issue. One 
drawback of NURBS is the lack of local refinement, thus making adaptive NURBS-based 
IGA impossible. To this end, T-splines have been recently developed to overcome the 
limitations of NURBS while making use of existing NURBS algorithms [5j. Another 
solution is the hierarchical refinement approach presented in [18]. These two options 
can be straightforwardly incorporated into our IGA-HMM. 

Table 1: Differences between NURBS-based Isogeometric analysis and Finite element analysis [TT] 





IGA 


FEA 


Geometry 


exact 


approximation 


Handling points 


control points 


nodal points 


Variables 


control variables 


nodal variables 


Basis 


NURBS 


Lagrange 


Basis interpolation property 


no 


yes 


Continuity 


easily controlled 


C°, fixed 


Refinement space 


hpk 


hp 


Positive basis property 


yes 


no 


Convex hull property 


yes 


no 


Variation dismishing property 


yes 


no 



Fig. [i] plots two meshes that consist of quadratic NURBS elements-one with C 1 
continuity across element edges (Fig. [4^i) and one with only C° continuity across element 
edges (Fig. [4Jd) . The latter was seamlessly obtained from the former by using the so- 
called knot insertion operation. Also clear from the referred figure is the fact that the 
number of dofs of the C 1 mesh is less than the number of dofs of the C° mesh. It is worthy 
noting that as a result of several decades of research, many efficient computer algorithms 
exist for the fast evaluation and refinement of NURBS [16] . Therefore, NURBS-based 
FEM is very efficient. 
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(a) C 1 quadratic NURBS mesh 



(b) C° quadratic NURBS mesh 



Figure 4: Quadratic NURBS meshes (2x2 elements): C 1 case (a) and C° case (b). Red 
dots denote the control points which generally do not locate inside the physical domain. 

3. An isogeometric analysis heterogeneous multiscale method (IGA-HMM) 

3. 1 . Model problems 

Let Vt be a domain in M. d (d — 1, 2, 3) with Lipschitz boundary dQ on which we impose 
the Dirichlet boundary conditions. Given /gL 2 (Q) as a source term, we consider the 
following classical highly oscillating coefficient second-order elliptic problem 

f-V-(a £ (x)W(x)) = /(x) infi, (] 
\ u e (x) = 0, on dn, { ) 

where 

• conductivity tensor a £ G (L°° (^)) dxci , and is uniformly elliptic and bounded, 

• e is a parameter which represents a fine scale characterizing the multiscale nature 
of a £ . 

and without loss of generality, we have omitted the Neumann boundary conditions. 

Equation ^ is popular in many practical problems such as solid mechanics, thermal 
conduction, and electrostatics, etc, where u can be the displacement field or electric 
field, etc, respectively, see Tab. [2j 
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Table 2: Significance of u and / in application [9 



Application Problem a £ 



f 



Solid mechanics 
Heat conduction 



Material tensor 
Thermal conductivity 
Acoustic conductivity 
Flow conductivity 
Flow conductivity 
Electrical conductivity 
Magnetical conductivity 



Temperature 
Displacement potential 
Pressure 
Velocity 

Electric potential 
Magnetic potential 




Mechanical force 
Heat flux 



Acoustic fluid 
Potential flows 
General flows 
Electrostatics 
Magnetostatics 



Particle velocity 
Particle velocity 
Fluxes 



Charge density 
Magnetic intensity 



In this paper, we confine ourselves to two dimensional heat elliptic problems, where 
u £ is the temperature field. We emphasize that the applications for other problems are 
similar and extension to three dimensions is straightforward [4]. 

From homogenization theory, we have known that the solution u £ of (1) (weakly) 
converges to u°, which is the solution of a so-called homogenization problem 



where the so-called homogenized tensor a is usually not available nor explicitly analyt- 
ically computed. 

The FE-HMM or IGA-HMM aims at finding the upscale solution u° without com- 
puting a explicitly. 

3.2. Drawbacks of the FE-HMM method 

One of the shortcomings of the FE-HMM is degree elevation. Due to the connection 
of different types of corner, center, or internal nodes, the implementation for degree 
elevation of standard FEM basis functions is, albeit possible, not an easy task. So far, 
to the best of the authors' knowledge, in the FE-HMM numerical literature, the highest 
order basis function has been used is cubic [17J. Another drawback of FE-HMM is that 
geometry errors will appear when applied to curved boundary domains. This is because 
in FEA meshing, polygons are used to approximately represent curved boundaries. To 
reduce this geometry error, we have to rely on sufficiently fine meshes, but this will lead 
to a high number of micro problems. These issues can be seamlessly solved when we 
use a new approach~the so-called isogeometric analysis heterogeneous multiscale method 
(IGA-HMM). Moreover, the IGA-HMM is able to utilize fc-refinement - higher order and 
higher continuity, which is a unique advantage of IG A [8] . 

3.3. The isogeometric analysis heterogeneous multiscale method (IGA-HMM) 

The crucial idea of the IGA-HMM is to replace the standard FEM basis functions 
used in FE-HMM by the NURBS basis ones, in both macro and micro scales. Before 



V • (a°(x)Vu°(x)) 

U°(x) 



0, x G oa 



(10) 
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describing the IGA-HMM method, we introduce the so-called macro and micro finite 
"patch" spaces in which we will work. 

Let ft be an open subset of R d , d — 1, 2, 3, and denote Y = (0, l) d the unit cube in 
M. d . We list here some familiar notations: 

1/2 

L 2 {Q)= {v.VL^^L: I \v\ 2 dx < oc \ , \\v\\ L 2 (n) = ( / \v\ 2 d* 



H 1 (fi) = {veL 2 (Q) : Vv e (£ 2 (ft)) d } , ||v||^i (n) = ^ 0| 2 + \Vv\ 2 )d^j 



1/2 



(fi) = {veH 1 (n):v = OondQ} ; 
^er(^) = {v e H X (Y) : v is Y- periodi^ . 



Remark 3.1. In this work, we assume NURBS with equal orders in two directions 
(p — q)- At the macroscale, NURBS of order p and at the microscale, NURBS of order 
q are utilized. This facilitates the derivation of error estimations presented later. 




Figure 5: IGA-HMM: from parameter mesh to physical mesh at the macro scale. For 
every Gauss point of the macro element K ) a micro domain Ks t is considered. The 
triangulation of K$ t is Th- Solution of the micro problem using the micro mesh Th yields 
the missing information at the macro Gauss point ■ 



3.3.1. Macro finite patch space 

We denote by 5 = ^ n+p+ i} , % = {771, rj m+p+ i} two knot vectors which 

define the index space, and F is a global geometry map function which parameterizes 
the physical space Q 

F : Q Q, 



''"For Y = (0, a function / : R d — >• R is said to be Y — periodic if and only if /(x + e^) = 
/(x) for all x G R d ,j G where {ej}j =1 is the canonical basis of R d . Thus we may view a 

function v G H^ er (Y) as a function from Y to R which belongs to iJ 1 (Y") and has the same trace on 
opposite faces of Y. 
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F is defined by the first equation of Eq. [8j The knot vectors 5 and T~L define, on the 
parameter space fi, a macro parameter mesh Th- Through the geometry map F, the 
parameter mesh Th become the physical mesh Th on the domain Vt. We denote by Hk 
the diameter of elements K e 7#, and set H := maxx^ff^. We refer to Fig. [5] for a 
graphics illustration. 

The macro finite patch space is then defined as 

SI (Q, T H ) = Hi (O) n spanji^. o F" 1 } . ., (11) 

where {Rij} are the two dimensional NURBS basis functions of order p defined on the 
parameter space Q. 

On each element K G Th, we consider {(jJk ^Ki)i = i ( p +i)2 the Gauss quadrature 
weights and integration point^Jas in the traditional FEM. In words, for a macro element 
of order we employ (p + 1) x (p + 1) Gauss points for numerical integration purposes. 
Optimal quadrature rules reported in [12] use less integration points than the standard 
quadrature rule (and hence results in less micro problems) are not used in this work. 

We set x.K t := F (x^) and consider in Q the following sampling domains 

K Sl =x Kl +6I (12) 
where / = (— |, with 5 > e. In words, the micro domain is a square centered on 

3.3.2. Micro finite patch space 

On each sampling domain K$ n we consider a micro triangulation Th induced by a 
geometry map which define K§ t 

F hA :K Sl ^K Sl , (13) 

where K$ l denotes the micro parameter domain which is (0, l) 2 . The partition Th consist 
of elements T of diameter We denote h := max^ G ^/i T . 

Now, the micro finite patch space is defined as 

S" (K Sl ,T h ) = W (K Sl ) n span{R htij o F^" 1 } . ., (14) 

where 



•'■In 2D, these Gauss points are defined as follows. First, they are defined in the usual parent domain 
[—1,-1] x [1,1]. Next, they are transformed to K. Note that this step does not exist in a standard 
FEM. Refer to [8 for details on integration issues of IGA. 
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• the Sobolev space W(Ks t ) depends on which kind of coupling (periodic or Dirichlet) 
between the macro and micro space functions is chosen [I]: 

W(K Sl ) = W l per {K 5l ) = i^v e Hl r {K Sl ) : jf vd* = J , (15) 
for periodic coupling, or 

W{K 5l ) = Hl{K 5l ) (16) 

for Dirichlet coupling; 

• {Rh,ij}ij are the two dimensional NURBS basis functions of order q defined on the 
sampling domain K$ t • 

Next, we describe the Isogeometric Analysis Heterogeneous Multiscale Method. 



3.3.3. The IGA-HMM method 

The goal of the IGA-HMM method is to give a numerical solution u H of the homoge- 
nization problem (10). The variational problem is presented after introducing the macro 
bilinear form and the micro problems. To facilitate the readers, the notations in this 
section follow reference pp. 



Macro-bilinear form We define the macro bilinear form Bh as 

B H (v H : w H ) := E Wl I ^( X ) V < ' V < rfx > ( 17 ) 
KeT H 1=1 1 Stl jK h 

where C is the number of quadrature points; v% and w\ § are the solutions of the 
following so-called micro-problems: 



Micro-problems Find v\ 8 (resp. w\ 6 ) such that 

'$u,k) eS*(K 5l ,%), 



and 



/ fl £ (x)v4 { • z h dx = 0, Vz h e S q (K 5l , T h ) 



(18) 



(19) 



where S q (Ks n 7h) is the micro finite element space defined in (14) and vf inK is the 
linearization of at the Gauss quadrature point which reads 



v Un,K, = V K t ( X xJ + VV H (x Kl ) • (X - X K J 



(20) 
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Variational problem The solution u H of the IGA-HMM is defined by the follow- 
ing variational problem: find u H G Sq (Q,Th) suc h that 

B H (u H , v H ) = (/, v H ) = / fv H dx, Vv H e S p (0, 7i) - (21) 
More details on the implementation can be found in [4] . 

3.4- A priori error estimates 

A priori error estimates for the NURBS-based isogeometric approach have already 
derived in [6], in which the authors showed that using NURBS basis functions of degree 
p in IGA, the error convergence rates are the same as those in the standard FEM with 
a polynomial of the same order p. This fascinating result allows us to predict that the 
a priori errors in the IGA-HMM follows the results derived from the FE-HMM [U [22] 
as shown in what follows. 



Let u°,u H be the solution of problems (10) and (21), respectively. Under some 



assumptions of regularity and periodicity on a £ ([lj, Proposition 14), using periodic 



coupling (15) and choose the size S of sampling domain as an integer multiple of e ) we 



have the following estimates for the solution obtained by the IGA-HMM method: 

h°-^ mm ^C^+(^y + ey (22) 

F-«X ( n ) <C'(^ 1 +(7) 29 + ^. (23) 

where the constant C depends on the shape of the domain Q and the shape regularity 
of the mesh, but does not depend on e ) the analytical solution u°, and the macro/micro 
mesh size. 

When the size s of the period is unknown, there is a so called cell resonance error 
due to a mismatch of the sampling domain size (e.g., 5/e ^ N) and to artificial boundary 
conditions arise, see [221 [23] . 

3. 5. Optimal refinement strategy for macro and micro meshes 

We denote by N maC) N mic respectively the macro and micro element per dimension 
of Q and the sampling domains. The corresponding macro diameter and scaled micro 
diameter are given by H = l/N maC) h = l/N miC) respectively. 



From estimates (22) and (|23j) we have the following optimal micro refinement strate- 
gies 

p+l r, 

Nmic = N mac ^ (optimal L strategy), (24) 



§e is very small, e.g. e « 10 6 , which makes the tensor oscillate very fast, but the u H given by the 
HMM method can capture the homogenized solution and "average out" these fluctuations. Also, for 
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p -1 

Nmic = N ma( M (optimal H strategy). (25) 



3.6. Total computational cost of the IGA-HMM and the choice of macro /micro degree 

Let N mac denote the number of macro elements per dimension, p, q the degree of 
macro and micro basis functions, respectively. Using the L 2 micro refinement strategy, 
we have: Total computational cost = number of macro elements x number of micro 
problem per macro element x number of dofs per micro problem. Thus, the total 
computational cost is 



Kac x (P + !) x (qN, 



P+l 
2q 

mac 



if 



2+ 



P+i 



From Eq. (26) we see that if the degree of the micro space is fixed at q 



(26) 

1, the total 



cost of the IGA-HMM will increase with the degree elevation of p: 



The total cost - 0(N££) (q = 1). (27) 

For p = 1 (this means linear basis functions are used in both macro and micro space), the 
complexity is of 0(N^ iac ) and then raise for each higher p, which is very inefficient. The 
total time to compute the solution will be very large and can be prohibited. Therefore, 
instead of fixing q = 1 at micro level, we choose q > p + 1. Then, 



The total cost - 0(N^ ac ). (28) 

In this case, the complexity of the IGA-HMM will be hold at 0(N^ ac ), for every p, 
which is acceptable and IGA-HMM computations can thus be run on normal personal 
computers or laptops. 



Remark 3.2. The formula (26) also holds in FE-HMM, but in this standard approach, 
it is very hard to increase the basis function degree when doing implementation due to 
the different connection between various kinds of nodes such as center, internal, edges, . . . 



4. Numerical examples 

In this section, we test the performance of the proposed IGA-HMM through three 
thermal problems. First, a problem on a square domain with analytical solution is 
analyzed with a linear basis for the micro element space. Convergence studies are per- 
formed for various macro NURBS basis functions orders (up to order five). Secondly, we 
solve a problem on a square domain that does not have an analytical solution. As the 
third problem, we consider a problem defined on a curved boundary domain geometry of 



tensor a £ (x) = a(x,x/e) with explicit scale separation, by collocating slow variable at the Gauss points 
and choosing 5 as an integer multiple of e, then the e part in the right hand side of Eqs. (22) and (23) 
will be "canceled out" HI. 
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which can only be exactly represented by high order NURBS. For the last two examples, 
high order (up to order five) NURBS elements are employed at both macro and micro 
scale which has not been done before in the framework of FE-HMM. To the best of our 
knowledge, for standard FEM basis functions, in the literature, the highest degree ever 
used in macro space is cubic, and in micro space is quadratic [TT] . 

The computations are performed on a desktop computer with CPU Intel® core i5 
2.8GHz processor. The present method has been coded in Matlab® language and for 
comparison purpose, the linear FE-HMM approach has also been implemented in our 
package. Due to its superior performance than the Dirichlet coupling, periodic coupling 
is chosen for all examples and 5 = e is used. A value of e = 10 -6 is used for all 
examples. Note that the IGA-FEM or FE-HMM captures the effective solution and 
is thus independent of e [4]. C° NURBS are used in the micro spaces unless otherwise 
stated. We recall that C° NURBS are similar to high order finite elements (C° continuity 
across element boundaries), refer to Fig. [31 



J^.l. IGA-HMM on domain with straight boundary 

We consider the following two-scale problem pQ 



(a(*) Vu e (x)) =1 in ft = (0,1)' 



(29) 



-V 

u e (x) = on dtt D := {x 1 = 0} U {x 2 = 1} 
[ n • (a Vu e (x)) = on dCl N := dn\dCl D , 

where a (*) = (cos (27r^) + 2)I 2 , x = (xi, x 2 ), and I 2 is the identity matrix. 

The corresponding homogenized tensor and solution can be computed analytically 



u°(x) 



X a 



2^3 



1 _dt_ 

k{t) 










'V3 0' 


2 




2 . 



where k(t) = cos27ri + 2. The analytical solution u is plotted in Fig. [6j 

We test the performance of IGA-HMM by considering the L 2 error and H 1 (energy) 



errors. 



In this example, we choose the degree q of the micro basis functions equal to 1. Hence, 
the micro element becomes as in FE-HMM. We study the convergence of the error 
measured in H 1 and L 2 norm with different NURBS basis order p in the macroscopic 



solver. Set q = 1 from (24), (25) we have 



^ (L 2 strategy), 



Nmic = N ma J (H 1 strategy). 
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Figure 6: Homogenized solution of the problem |4.1 



4.1.1. IGA-HMM: convergence rate test 

First, we study the L 2 micro refinement strategy. From Fig. u\ we see that for linear, 



quadratic, cubic, and quart ic NURBS, the convergence rates of the errors in L 2 norm 
when using the L 2 micro refinement strategy have the slopes of 2,3,4,5, respectively. 
These match well the predicted theory. The details of the errors are given in Tab. [3j 



Mesh 


Method 


2x2 


4x4 


8x8 


FE-HMM a 


0.28950 


0.08360 


0.02100 


p=2 (IGA-HMM) 


0.06092 


8.36e-3 


1.03e-3 


p=3 (IGA-HMM) 


0.03452 


2.13e-3 


1.34e-4 


p=4 (IGA-HMM) 


0.01443 


5.34e-4 


1.66E-05 



In this work, the standard QA elements are used for FE-HMM in both macro and 
micro finite element space. 



Table 3: L 2 error of problem 4.1 , using the L 2 micro refinement strategy. 
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slope=3.02 



slope=3.99 
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10' 
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Figure 7: L error of problem 4.1, using the L micro refinement strategy 



For the H 1 error when using the L 2 micro refinement strategy, from Fig. |8j we see 
that the convergence rate is slightly higher than the theoretical value. This is due to 
the problem dependence. In the later examples, this phenomenon does not occur. We 
list in Tab. H]the detailed errors. 



Method 
FE-HMM 
p=2 (IGA-HMM) 
p=3 (IGA-HMM) 
p=4 (IGA-HMM) 



2x2 
0.50559 
0.06092 
0.03435 
0.01443 



Mesh 
4x4 
0.25168 
8.36e-3 
2.13e-3 
5.34e-4 



8X8 
0.12527 
1.03e-3 
1.34e-4 
1.66E-05 



Table 4: H 1 error of problem 4.1 using the L 2 micro refinement strategy. 




Figure 8: H 1 error of the problem AA_, using the L 2 micro refinement strategy 
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We next continue to test the problem for the H 1 micro refinement strategy. Fig. [9] 
and Tab. [5] demonstrate that for linear, quadratic, cubic, quartic, and quintic NURBS, 
the convergence rates of errors in H 1 norm are 1, 2, 3, 4, 5 and in L 2 norm (see Fig. 10 
and Tab. [6]) are 1,2,3, 4, 5, respectively. It is clear that the results match the theoretical 
ones. 



Mesh 



Method 


2x2 


4x4 


8x8 


16x16 


32x32 


64x64 


FE-HMM 


0.50559 


0.27581 


0.13888 


0.06725 


0.03443 


0.01770 


p=2 (IGA-HMM) 


0.07647 


0.02964 


8.31e-3 


2.13e-3 


5.34e-4 




p=3 (IGA-HMM) 


0.06439 


8.34e-3 


1.03e-3 


1.34e-4 






p=4 (IGA-HMM) 


0.03108 


2.13e-3 


1.34e-4 








p=5 (IGA-HMM) 


0.01414 


5.34e-4 


1.66e-05 









Table 5: H 1 error of the problem 4.1 , using the H 1 micro refinement strategy. 




Figure 9: H 1 error of problem 4.1, using the H 1 micro refinement strategy. 



Mesh 



Method 


2x 2 


4x4 


8x8 


16x16 


32x32 


64x64 


FE-HMM 


0.28950 


0.16737 


0.07342 


0.02809 


0.015260 


8.51e-3 


p=2 (IGA-HMM) 


0.07646 


0.02968 


8.31e-3 


2.13e-3 


5.34e-4 




p=3 (IGA-HMM) 


0.06471 


8.35e-3 


1.03e-3 


0.00013 






p=4 (IGA-HMM) 


0.03120 


2.13e-3 


1.346-4 








p=5 (IGA-HMM) 


0.01413 


5.34e-4 


1.66e-05 









Table 6: L 2 error of the thermal square problem, using the H 1 micro refinement strategy. 
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slope=1 .19 



slope=1 .84 
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Figure 10: L 2 error of problem 4.1, using the H 1 micro refinement strategy. 



4.1.2. IGA-HMM vs FE-HMM: performance comparison 

In this section we compare the CPU-time of the IGA-HMM and that of the FE-HMM 
needed to obtain solutions with various accuracies: 10 -2 , 2x 10 -3 , 5x 10 -4 , 10 _4 ,and 10 -5 
in H 1 norm of the error, using the H 1 micro refinement strategy. The results are pre- 
sented in Tab. [7j We can see that for the same accuracy, the IGA-HMM outperforms the 
(linear) FE-HMM in terms of CPU-time needed to compute the solution. For example, 
to obtain the H 1 error less than 10 -2 , the FE-HMM needs about 1300 second, while the 
IGA-HMM needs only approximately 5 seconds and 2.5 seconds for quadratic and cubic 
NURBS, respectively. Similarly, to obtain the error in H 1 norm approximately 2 x 10 -3 , 
the FE-HMM needs more than 80000 seconds, while the IGA-HMM needs only approx- 
imately 83 seconds and 16 seconds for quadratic and quart ic NURBS, respectively. We 
can see that the CPU time of the IGA-HMM is significantly reduced in comparison 
with the FE-HMM thanks to the high order NURBS basis functions employed at the 
macroscale which allows very coarse macro meshes to be used. It should be noted that 
our Matlab implementation was not optimized and hence performance can be further 
improved. 



4.2. A high- order approximation of IGA-HMM in both macro and micro 
levels patch spaces 

In the previous example, we have already seen that the IGA-HMM performs much 
better than FE-HMM when using high order NURBS (p > 2) in macro patch space 
and NURBS with degree q = 1 in micro space. In this section, we continue to test the 
IGA-HMM with high order NURBS in both micro and macro space. To this end, we 
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method 



mesh 



H 1 err 



CPU(s) 



T — \~i — 1 TTH T~\ T 

FE-HMM 


128 x 128 


8.67e-3 


1359 


O / T A TT1\ /TTV /T\ 

p=2 (IGA-HMM) 


8x8 


8.31e-3 


4.9 


/ T" ✓ — ^ A ~TT~A IV "TV /T \ 

p=3 (IGA-HMM) 


4x4 


8.35e-3 


2.5 


"I ^ T — ^ TT1\ /T1\ T 

FE-HMM 


512 x 512 


2.17e-3 


>80000 


o /t/^i a tti\ /ri\ /r\ 

p=2 (IGA-HMM) 


16 x 16 


2.13e-3 


83 


A / T /'"N A TTT\ /TT\ T \ 

p=4 (IGA-HMM) 


4x4 


2.13e-3 


16 


FE-HMM 


2D48 x 2048 






p=2 (IGA-HMM) 


32 x 32 


5.34e-4 


1433 


p=5 (IGA-HMM) 


4x4 


5.34e-4 


102 


FE-HMM 


8192 x 8192 


1.35e-4 


>3e8 (est) 


p=3 (IGA-HMM) 


16 x 16 


1.34e-4 


3064 


p=6 (IGA-HMM) 


4x4 


1.34e-4 


673 


FE-HMM 


65536 x 65536 


1.69e-5 


>2el0(est) 


p=5 (IGA-HMM) 


8x8 


1.66e-5 


6.06e4 



Table 7: Performance comparison between the FE-HMM and the IGA-HMM (both using linear basis 
functions in micro space). 



consider the following two-scale problem: 

-V • (a (x, *) V« e (x)) = 1 in ft = (0, l) 2 , 
< u £ (x) = on dft D := {x l = 0} U {x 2 = 1} , (30) 
n • (a (x, *) V« e (x)) = on dtt N := dQ\dQ D , 
where the conductivity tensor is given by [3]: 



a (x, — 



x\ + 0.2 + (x 2 + 1) (sin (2tt^) + 2) 







The corresponding homogenized tensor is given by 

" ( f 1 dm y 1 

a °(x) = ^ x l+ - 2 +( a; 2+l)(sm(27ry 1 )+2) 



0.05 + {x x x 2 + 1) (sin (2tt^) + 2) 

(31) 









dm 



i|+0.05+(xix 2 +l)(sin(27ry2)+2) 



(32) 



This problem does not have an explicit analytical solution, thus for a reference solution, 
we solve the homogenized problem with the homogenized tensor a on a fine mesh of 
500 x 500 quintic elements using the standard NURBS-based finite element method, see 
Fig. 11 The relative errors are computed using the following formula: 



err H i 



\Uh u ref\\ ffi(Q) 



\u. 



re/ llifi(n) 



err L 2 



\Uh U ref\\ L 2^ 



U, 



"refWn 1 ^) 
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For the macro space we use NURBS of degree ranging from 1 to 3, and for the micro 



space, from the choice q > p + 1 as discussed in Eq. (28), we use the NURBS function 
of the degree fixed at 5 (q should not be too large, because when the micro mesh is too 
coarse, the results will not stable). The results are presented in Tabs. M and [9} 




o o 



Figure 11: Homogenized solution of the problem 4.2 









Mesh 






Degree 


4x4 


8x8 


16 x 16 


32 x 32 


64 x 64 


FE-HMM 

p=2 

p=3 


0.24884 
6.97 e-3 
6.73 e-4 


0.12443 
1.71 e-3 
9.01 e-5 


0.06223 
4.26 e-4 
1.20 e-5 


0.03112 
1.06 e-4 
1.59 e-6 


0.01557 
2.66e-5 
2.07e-7 




Table 8: H 1 


error of problem 4.2, using the L 2 


micro refinement strate 


gy 








Mesh 






Degree 


4x4 


8x8 


16 x 16 


32 x 32 


64 x 64 


FE-HMM 

p=2 

p=3 


0.06063 
7.66 e-4 
6.51 e-5 


0.01514 
8.95 e-5 
4.10 e-6 


0.00379 
1.09 e-5 
2.70e-7 


0.00095 
1.37 e-6 
3.24 e-8 


0.00024 
1.72e-7 
l.lle-9 



Table 9: L 2 error of problem 4.2, using the L 2 micro refinement strategy. 
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Figure 12: H 1 error of problem 42, quintic NURBS in micro space and L 2 refinement 
strategy are used. 



10" 
10" : 

io" ; 

10"' 
10"' 
10"' 
10" 
10"' 
10" ! 









^^^^ slope=2 




^^K. slope=3 






— » — deg=1 


slope=4 


► deg=2 




— ★ — deg=3 





10 



10' 

1 / Umax 



10 



Figure 13: L 2 error of problem 42, quintic NURBS in micro space and L 2 refinement 
strategy are used. 



Figures [12] and [13] shows that the convergence rate of the IG A-HMM is optimal and 
matches the theory suggested by FE-HMM: 1,2,3 in H 1 error and 2, 3, 4 in L 2 error for 
linear, quadratic and cubic NURBS in the macro space, respectively. From Tabs j8] and [9j 
we can see that the IGA-HMM is able to obtain a very high accuracy, up to 2.07e — 7 and 
1. lie — 9 in H 1 error and L 2 error, respectively. Such accuracy is almost prohibited in 
FE-HMM, as pointed out in the previous example. Also in comparison with the results 
obtained in the previous example, we can see that the high order approach in both micro 
and macro spaces is more effective than when we only use high order in macro and linear 
basis functions in micro spaces, which is only up to 10 -6 in L 2 error. 
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4.3. IGA-HMM applied for a curved boundary domain 



In this problem, we test the convergence of the IGA-HMM method on a curved 
boundary domain, where traditional FEM methods cannot exactly represent the bound- 
ary. Let be the quarter circular annulus of internal radius r = 1 and external radius 



R — 2, which lies in the first quadrant of the Cartesian plane, see Fig. [15} Homogeneous 
Dirichlet boundary condition is imposed on the whole boundary. Let the heat source be 



/ = 



2x^x 2 (18V£i 2 + £2 2 + 4V3-24) + 2x x x^ (9 V3 ^x x 2 + £ 2 2 - 12 ^3 ■ 



(xi 2 + x 2 2 ) 



2\ 3 



We consider the following problem: 



-V • (a V« £ (x)) = / in ft, 
m £ (x) = on 9ft 

where a (|) = (cos (27r^-) + 2)I 2 , x = (xi, x 2 ), and I 2 is the identity matrix. 



(33) 



(34) 



The analytical homogenized solution u is given by (see also Fig. 14) 



u 



X\X2 (4 — 6\/^l 2 + ^2 2 ) 



(35) 



X\* + x 2 * 

First, we consider the case with linear basis function at micro space and for three macro 



meshes as shown in Fig. 16, The geometry data for the coarsest mesh (single quadratic 



element) is given in Tab. 10, Other meshes of different densities and basis orders 

PI 



are constructed seamlessly from this coarsest mesh. Tables [TTJ [12] and Fig. [171 
demonstrate that the convergence rate of the error is optimal and matches the theory 




Figure 15: Coarsest mesh (one single element) of the domain of problem 4.3 and its 
control net. 




i 


1 


2 3 


4 


5 6 


7 


8 9 


Pi 


(0,1) 


(1,1) (1,0) 


(0,1.5) 


(1.5,1.5) (1.5,0) 


(0,2) 


(2,2) (2,0) 


Wi 


1 


l/y/2 1 


1 


l/y/2 1 


1 


l/y/2 1 



Table 10: Control points and weights for the coarsest mesh of the computational domain of problem 



4.3 The corresponding knot vectors are S = {0, 0, 0, 1,1,1},% = {0, 0, 0, 1, 1, 1}. 



Degree 




Mesh 




2x2 


4x4 


8x8 


p=2 


0.06920 


0.01682 


0.00439 


p=3 


0.02575 


0.00423 


0.00052 


p=4 


0.01567 


0.00107 


6.60E-05 


p=5 


0.00713 


0.00026 


8.14E-06 



Table 11: H 1 error of the thermal quarter annulus problem, using linear basis function in micro space 
and the H 1 micro refinement strategy. 
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Figure 17: H 1 error of the thermal quarter annulus problem, using linear basis function 
in micro space and the H 1 micro refinement strategy. 



Degree 




Mesh 




2x2 


4x4 


8x8 


p=2 


0.05895 


0.01509 


0.00400 


p=3 


0.02501 


0.00401 


0.00049 


p=4 


0.01514 


0.00102 


6.39E-05 


p=5 


0.00699 


0.00026 


7.91E-06 



Table 12: L 2 error of the thermal quarter annulus problem, using linear basis function in micro space 
and the H 1 micro refinement strategy. 



We emphasize that in FE-HMM approach, it is very hard to use high order elements 
for problems with curved boundary due to the requirement of the mesh quality: the 
partition of the domain must be fine enough so that the curved boundary can be ap- 
proximately good enough. As a consequence, the computational cost will be extremely 
high, if not prohibited. 

Next, we use high order NURBS basis functions in both micro and macro spaces. 
Macro space utilizes NURBS of degree ranging from 2 to 5, while in micro space quintic 
NURBS is used. Here, we follow the L 2 micro refinement strategy. The results in Tabs. 



13, 14, and Figs. 19, 20 show that this approach greatly outperforms the standard 



FE-HMM in terms of convergence rate, accuracy as well as efficiency. In Tab. [15] and 
Fig. [2TJ the results also confirm the nearly optimal computational cost given by Eq. [26| 

It is worthy noting that for problems with curved boundaries as this example, it is 
very hard for FE-HMM to obtain accuracy of magnitudes 10 -6 , 10 -7 , 10 -8 in H 1 error 
or 10 -7 , 10 -8 , 10 -9 in L 2 error, because of the high number of macro elements required 
to approximate the boundary. Furthermore, elevating the degree of basis functions in 
standard FEM is, albeit possible, not a flexible task. Here, in IGA-HMM with high 
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Figure 18: L 2 error of the thermal quarter annulus problem, using linear basis function 
in micro space and the H 1 micro refinement strategy. 

order NURBS basis functions in both macro and micro space we can obtain these high 
accuracies very easily: about 1000 seconds with mesh 8 x 8 (p = 5) for the accuracy of 
10 -6 in H 1 error and 2 x 10 -7 in L 2 error; for accuracy of 2.7 x 10 -8 in H 1 error and 
2.7 x 10 -9 , it takes about 9000 seconds with a mesh of 16 x 16 (p = 5). 



Mesh 



Degree 


2x2 


4x4 


8x8 


16 x 16 


p=2 


0.04501 


0.00724 


0.00164 


0.00040 


p=3 


0.00450 


0.00116 


1.13E-04 


1.33E-05 


p=4 


0.00241 


2.72E-04 


1.04E-05 


5.79E-07 


p=5 


0.00012 


3.80E-05 


1.11E-06 


2.71E-08 


Table 13: H 1 


error of the thermal quarter annulus problem, usin;; 


5 NURBS of dej 


$ree 5 in micro space 






Mesh 






Degree 


2x2 


4x4 


8x8 


16 x 16 


p=2 


4.27 


17.01 


68.02 


512.62 


p=3 


7.96 


32.33 


235.54 


1638.44 


p=4 


12.48 


49.84 


374.17 


2592.46 


p=5 


18.84 


141.07 


956.76 


9055.52 



Table 15: CPU time (s) solving the thermal quarter annulus problem, using quint ic NURBS in micro 
space. 
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Mesh 



Degree 


2x2 


4x4 


8x8 


16 x 16 


p=2 


0.02784 


0.002084 


0.00022 


2.58E-05 


p=3 


0.00211 


0.00041 


1.89E-05 


1.08E-06 


p=4 


0.00108 


1.10E-04 


1.84E-06 


5.78E-08 


p=5 


4.43E-05 


1.19E-05 


2.08E-07 


2.71E-09 


Table 14: L 2 error 


of the thermal quarter 


annulus problem, using 


• NURBS of dej 


;ree 5 in micro space. 
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Figure 21: CPU time (s) solving thermal quarter annulus problem, using quintic NURBS 
in micro space. 



Hitherto, C° NURBS are used in the micro spaces. In what follows, we test the 
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Figure 20: L 2 error of the thermal quarter annulus problem, using quintic NURBS in 
micro space. 



performance of the IGA-HMM when utilizing high order NURBS with C p ~ l continuity 
in the macro space and with C q ~ x continuity in the micro space. The errors in H 1 norm 
and L 2 norm are given in Tabs. 16 and 17, respectively. In this case, the performance 
of the method is worst than the performance when C° NURBS is used in micro spaces, 
and the convergence rate is not optimal, see Fig. [22} The reason for this phenomenon, 
to our knowledge, is that in the micro space the material is strongly heterogeneous as 
reflected in a (*). As a consequence, high order C q ~ l continuity NURBS basis functions 
does not capture well these heterogeneities as the C° ones, although the accuracy is still 
acceptable in comparison with that obtained with the standard FE-HMM. 







Mesh 






Degree 


4x4 


8x8 16x16 


32x32 


64x64 


p=2 
p=3 


0.045027 
0.004546 


0.007264 0.001738 
0.001299 2.7E-4 


4.67E-04 
1.0E-04 


2.70E-4 
1.0E-4 



Table 16: H 1 error of the thermal quarter annulus problem, using quintic NURBS in micro space with 
C° continuity 



Mesh 

Degree 4x4^ 8x8 16x16 32x32 64x64 

p=2 0.027884 0.002179 0.000553 0.000221 0.000230 

p=3 0.002213 0.00070 0.000231 9.27E-05 9.42E-5 



Table 17: I? error of the thermal quarter annulus problem, using quintic NURBS in micro space with 
C° continuity 



28 



10° 10 1 10 2 10° 10 1 10 2 

l/H max l/H max 

Figure 22: H 1 and L 2 errors of the thermal quarter annulus problem, using quintic 
NURBS in the micro space with C p ~ l continuity. 

5. Conclusions 

We have, for the first time, presented an efficient isogeometric analysis heterogeneous 
multiscale method (IGA-HMM) for elliptic homogenization problems. The method is 
capable of capturing the exact geometric representation and is very flexible regarding to 
refinement and degree elevation by using the NURBS basis functions in both macro and 
micro levels. As a result, the high-order IGA-HMM macroscopic and microscopic solvers 
are designed simply and effectively. A priori error estimates of the discretization errors 
and optimal micro refinement strategies were derived. The provided numerical results 
showed that the IGA-HMM achieves high order accuracy with optimally convergence 
rate for both L 2 and H 1 micro refinement strategies as in the theory of the standard 
FE-HMM. We have also observed that C°, not C p ~ l , high order NURBS should be 
used at the microscale although this point deserves a further study. The importance 
is that this approach is very similar to the standard FE-HMM, this means that the 
implementation is very simple, the coding framework of FE-HMM in [4J can be re-used, 
while the accuracy and the flexibility are remarkably improved. These advantages make 
the IGA-HMM an alternative method to solve homogenization problems, beside some 
new methods have recently been published [TT1I3]. The method shows potential in solving 
high order homogenization problems due to the arbitrary smoothness of the NURBS. 

In order to apply the IGA-HMM for problems with a more higher complexity, a com- 
bination of IGA-HMM with techniques introduced in [15] can be considered. In this way, 
not only the number of macro elements will be reduced but also that of micro problems. 
As a result, the computational cost will dramatically be minimized. This will be our 
forthcoming development. 
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